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Segmentation and Restoration of Images on Surfaces by Parametric 
Active Contours with Topology Changes 
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Abstract 

In this article, a new method for segmentation and restoration of images on two-dimensional 
surfaces is given. Active contour models for image segmentation are extended to images on 
surfaces. The evolving curves on the surfaces are mathematically described using a parametric 
approach. For image restoration, a diffusion equation with Neumann boundary conditions is 
solved in a postprocessing step in the individual regions. Numerical schemes are presented 
which allow to efficiently compute segmentations and denoised versions of images on surfaces. 
Also topology changes of the evolving curves are detected and performed using a fast sub¬ 
routine. Finally, several experiments are presented where the developed methods are applied 
on different artificial and real images defined on different surfaces. 


Keywords: Image segmentation, images on surfaces, evolving curves on surfaces, active con¬ 
tours, parametric method, Mumford-Shah, Chan-Vese, topology changes, triple junctions, image 
restoration, finite element approximation. 


1 Introduction 


We study the problem of segmentation and restoration of images defined on two-dimensional 
surfaces. The Mumford-Shah functional |23] can be extended and reformulated for image data 
given on surfaces. Segmentation aims at dividing an image on a surface in characteristic regions, 
for example in regions of similar gray-value or similar color. The objective of image restoration is 
to denoise the original image while preserving sharp edges in the image. 

Active contours |l5j originally developed for planar images can also be used to segment images 
on surfaces. We let one or more time-dependent curves T(t), t E [0,T], evolve in time on a surface 
Met 3 according to a flow such that the curves are attracted by edges or region boundaries. 

Existing studies and works on evolution of curves on surfaces and on active contours for images 
on surfaces differ in the way how the surface and the curves are described mathematically. A 
geometric scale space for images on parametric surfaces is introduced in 16 using level sets. The 


image is handled implicitly by considering its iso-gray levels. In 28 , Spira and Kimmel consider 


flows of curves on parametric surfaces and perform edge detection using a variant of the geodesic 
active contours method |6 for images on surfaces. In |28], the authors restrict on surfaces which 
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have a global parameterization. Evolution equations are solved with the level set method [25] by 
considering the pre-image of the curve resulting in a planar curve in the parameterization plane. 
However, this approach has intrinsic disadvantages since the pre-image of the curve is used. In |17], 
drawbacks of this approach concerning the scaling behavior are discussed. Another drawback is 


the fact that the method does not allow to incorporate a balloon term 17 


In [10], a different method to represent the surface is proposed, where the surface is modeled by 
the zero level set of a three-dimensional function. The zero level set of a second three-dimensional 
function, which is time-dependent, is used additionally. The curve is represented by the intersection 
of the two zero level sets. This approach is used in 17 and [331 for image segmentation with 


geodesic active contours. A drawback of the method is that a one-dimensional curve evolution 
problem is extended to a three-dimensional problem. To reduce the effort, a new narrow band 
technique is given in 1171. 

In 1291, image segmentation is performed using the Chan-Vese |9| model for images on surfaces 
with level set methods in combination with a so-called closest point method. One iteration step 
of the Chan-Vese model in a small 3D neighborhood of the surface is computed followed by an 
interpolation step. 

Flows of the form V n — km + F are studied in 


22 , where V n is the velocity contribution 


normal to the curve, but tangential to the surface, km is the geodesic curvature of the curve, and 
F is an external forcing term. The authors restrict on graph surfaces and solve a system of partial 
differential equations for the parameterization x, the tangent angle, the geodesic curvature km 
and for \\x p \\. 

In this work, we make use of a concept developed in |3|, where Barrett et al. consider flows of 
curves on surfaces like geodesic curvature flow, geodesic surface diffusion and geodesic Willmore 
flow of curves. Here, we apply the methods of 131 to image segmentation applications resulting in a 
curvature driven flow: The normal velocity is the sum of the geodesic curvature km weighted with 
a parameter a and an external forcing term. The forcing term is designed for image segmentation 
based on an extension of the Chan-Vese method to images on surfaces. As in |3|, a parametric 
approach is used to describe the curve. The surface is given implicitly, as zero level set of a smooth 
function $. However, for practical image segmentation, the function $ need not be explicitly 
known. Our method requires only a normal vector to the surface at each point on the surface. The 
developed segmentation technique can handle multiple phases and networks of curves including 
triple junctions. 

The second image processing task we consider in this article is image restoration. To denoise a 
given image, we solve a diffusion equation on the surface AT For related works on surface partial 
differential equations, we refer to 1131 and the references therein. For a practical application 


see 


11 , where data of the Earth’s surface captured on board of satellites are filtered by nonlinear 


diffusion. An implicit representation of the surface is used in |5], where the surface is embedded 
as zero level set of a higher dimensional function and the partial differential equations are solved 
on a fixed Cartesian grid using a special embedding function. 

Methods to solve total variation problems on surfaces are proposed in 


18 . In detail, the 


method of |27] for denoising images on surfaces and the method of |7| (a convexified Chan-Vese 
model) for segmentation of images on surfaces are generalized. A direct, called intrinsic , approach 
is pursued: Lai and Chan |18j perform the resulting calculations directly on the given surface by 
using differential geometry techniques and finite elements for images given on triangulated surfaces. 
They also provide an overview and a comparison of several approaches for variational problems 
on surfaces (level set methods, parametric surfaces and direct/intrinsic methods). Total variation 
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based image restoration and segmentation are also considered in |32|, where the authors propose 
an extension of the augmented Lagrangian method (see e.g. |3l|) for scalar and vectorial total 
variation problems to images on surfaces. 

In this article, we will perform restoration of images on surfaces by considering an extension 
of the Mumford-Shah problem [23]. The image restoration is performed as a postprocessing step 
after the segmentation. We solve a diffusion equation with Neumann boundary conditions in the 
already segmented regions. Thus, the image is not smoothed out across the regions boundaries. 

For both segmentation and restoration we present efficient numerical schemes. Further, topol¬ 
ogy changes of the parametric curves are detected efficiently. In [4], we used and extended a method 
of |21 for efficient detection of topology changes, such that also topology changes involving triple 
junctions can be handled during the curve evolution. We will extend this idea to curves on surfaces 
to detect several topology changes including splitting and merging of curves, and creation of triple 
junctions. 

In practical applications, surfaces are often not given as smooth surfaces but in a discrete 
form, for example as triangulated surfaces. We present all necessary computational aspects when 
applying the segmentation and restoration models to real data. 

In summing up, we develop a novel scheme for both segmentation and restoration of images 
defined on surfaces. Using our parametric approach, the evolution of curves is a one-dimensional 
problem; the postprocessing image diffusion is a two-dimensional problem computed only once after 
the segmentation has been finished. Compared to other approaches in the literature, where the 
curve is embedded in a higher-dimensional space, our method is very efficient from a computational 
point of view. 

2 Segmentation and Restoration of Images on Surfaces 

2.1 Preliminaries 

Let M C K 3 be a smooth two-dimensional manifold. We assume that we can describe M by the 
zero-level set of a function <L: 

M = {zeR 3 : $(z) =0}, (1) 

where $ G C 2 (R 3 ,R) is a function with ||V$(z)|| > 0 for z G M. A unit normal vector field n<$> 
on M is given by n$(z) := V$(z)/||V$(z)|| for z G M. 

Let T C M be a curve on A4 parameterized by x : / —> A4, where / is a one-dimensional 
reference manifold, for example the unit interval / = [0,1] for open curves or R/Z for closed 
curves. We define : / —» M 3 such that z7$(p) := n$(x(p)) is the surface normal evaluated at 
x(p) for p G /. Further, we define Vm : I R 3 by i' m(p) %s(p) x z%(p), where s denotes 

the arc-length. We note that i7 m is perpendicular to x s , i.e. normal to the curve, but lies in the 
tangent space to M. Figure [l] illustrates a possible surface M, a curve T and the vector fields 
x s and I'm- 

Further, the vector x ss which is perpendicular to x s can be written as the sum of its component 
in z/^-direction and its component in zTy^-direction. This motivates the following definition |3|: 

Definition 1. We define the geodesic curvature km \ I R and the normal curvature : / —>► R 

by 

f^J\4 = %ss • FA/p %ss • (2) 
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Figure 1: Illustration of the vector fields x s , V<§> and Vm — %s x 

As a consequence, the vector field x ss can be expressed as 

%ss = KM?M + (3) 

2.2 Active Contours on Surfaces Based on Extensions of the Mumford-Shah 
and Chan-Vese Models to Images on Surfaces 

Let uq : Ai [0,1] be the intensity function of an image given on a surface Ai. For segmentation 
and restoration, we consider an extension of the Mumford-Shah functional [23] from the planar 
case to the case of images and curves on surfaces. For that, we consider the following minimization 
problem: 

Find a union of curves V = Ti U ... U rN c C Ai and a piecewise smooth approximation 
u : Ai M of the original image uq such that 

E MS (u,T) = a\r\+ [ \\V M u\\ 2 dA +X [ (u 0 -u) 2 dA (4) 

Jm\t Jm 

is minimized. Here, cr, A > 0 are weighting parameters, |T| denotes the total length of the curves, 
and V m u is th e surface gradient of ix, also called tangential gradient, cf. |13|. Further, d A denotes 
the area element and || . || the Euclidean norm. 

We first consider the case of one closed curve T on M without any self-intersection which is 
homotopic to a point. Then, the curve divides Ai in two disjoint regions and such that 

M = urufi 2 . 

The indices of the regions k — 1,2, are chosen such that um points from to 

For segmentation, we consider a piecewise constant approximation with u\Q k = G R, k = 1, 2. 
The Mumford-Shah functional for images on surfaces 0 reduces to 

E(T,a,C 2 ) = cx|r| + A / (uq - ci) 2 d^4 + A f (uq - C 2 ) 2 dA. (5) 

J Oi J O2 

Similar to this functional, we can consider the analogue of the planar Chan-Vese model |9| for 
images on surfaces: 


E(r,ci,c 2 ) 


a\T\ + fi / 1 d A + Ai 


fi dA + A 2 / / 2 dA, 
J O 2 


( 6 ) 
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where a, Ai, A 2 > 0, p > 0 are weighting parameters. Similar to the planar Chan-Vese method |9|, 
the function /&, k — 1,2, is defined by 


— OoO) Oc) i ^ ^ •A'L 


( 7 ) 


The model © with fi = 0 and Ai = A 2 is the piecewise constant case © of the Mumford-Shah 
model. 

Fixing now the curve r in Q, we obtain the following condition for the coefficients c&, k = 1,2: 


Ck 



( 8 ) 


Let x : / = M/Z —Ad be a smooth parameterization of T. 

We now derive a flow for image segmentation as gradient flow using methods of the theory of 
calculus of variations. Since we consider curves on the surface Ai, the variations are restricted to 
he on Ai. Therefore, as in 131, we consider the variations to be elements of 

V_<£ = [fj : / —> R 3 : ff is smooth and ff. — 0} (9) 

and define for functions 77 , y : / —)► M 3 the following inner product 

%X)2 ,A4,nor • — l PmV'PmX^s, (10) 


where Pm is the projection onto the part in direction 1 7^, i.e. PmP — (v- Fm)fm for ff: I M 3 . 

Fixing the coefficients ci,C 2 , we consider for ff G V_§ a variation y : I x (—eo,eo) Ai with 
y(p,0) = x(p) and y e (p, 0) = ff(p). 

Let PcMbe the image of y (., e), let Of, 0| C Ai be regions such that Ai = Of U T e U 0|. 
Further, let ^m(p) G T^ pe ^Ai be a vector in the tangent space to Ai in y(p,e) defined by 
1 'm(p) = Vs(pa) x n$(y(Pi 6 ))* The vector field V ^(p) points in the direction Of. The vector 
Vm(p) is normal to T e , but lies in the tangent space of Ai at the point y(p, e). 

We define 


(SE(F)M := ^ 


e=0 


f cr f 1 ds + /i f 1 d A + Ai f f\ d A + A 2 f $2 d-A j 

\ J r e Jni Jni Jn e 2 J 


on noting that y (., e) and thus T e , Of and 0| depend on ff. We compute 


mmv) = ^ 


(a[ \\y P \\dp + p [ lcL4 + Ai [ fi dA + X 2 [ / 2 d^4 ) 
e =o \ Ji Jn\ Jn\ 2| y 


= cr 


y P 


y pe dp + p J{-y e .V e M )\\y p \\dp+\i J^f 1 (y)(-y e .iJ e M )\\y p \\dp + 


e=0 


1 II v P 

+A 2 J h{y)(ye ■ "mWpW d P 
= a J x s .fj s ds + J (-p - A 1/1 + A 2 / 2 ) VM -yds 

= / (~ax ss + (-p - A 1/1 + A 2 / 2 )^x) -yds 


= J ~ P ~ A 1/1 + A 2 / 2 ) v M .rjds. 


( 11 ) 
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Here, we used a transport theorem for curves on surface s [14] . We applied integration by parts for 
the second last identity. The last identity follows from ^ and ff . — 0. 

A time-dependent function x : / x [0,T] -G Ai with x t ( . ,£) G V_^ is called a solution to the 
gradient flow equation if 

(xt, ff)2,M,nor = ~ (SE( f)) (ff) (12) 


holds for all ff G V_$. 

Let 77 G V_$. We conclude from © and ( [T2] ) on noting that Vj^\ . ff = i7 m • PmV 


PMXt = - (-CTKX - n - Al/l + A2/2) VM- 


(13) 


Let V n := xt • Am denote the velocity in direction 1 7^, also called normal velocity. Then Pjv[^t — 
V u Vm and consequently the equation above leads to 


V n = cr nj\/i + F, 


(14) 


where F is given by 

F(z) = fj, + Xifi(z) - X 2 f 2 (z) = fi + Xi(u Q (z) -ci ) 2 - X 2 (u 0 (z) - c 2 ) 2 . (15) 

We rewrite equation ( [l4| ) to a scheme for x : / x [0,T] -G M 3 and : I x [0,T] -G M. 

We assume that x(p, 0) lies on M. To force the curve to stay on the manifold Ad, the velocity in 
direction normal to the surface xt -n® must be zero (i.e. xt G V_$). We thus have the following 
scheme: 

Let x(I, 0) = T(0) C M. For t G (0, T], find x{ ., t) : I -G M 3 and k>m( • , t), /s$( ., £) : / Gl 
such that 


Xt-VM = ^AM + F, (16a) 

xt.v<§> — 0 , (16b) 

^55 = AM AM + (16c) 


2.3 Multiphase Image Segmentation with Possible Triple Junctions 

We extend the above presented two-phase segmentation with a single closed curve to more general 
situations. We consider a curve network with closed and open curves which partition the image 
domain in Nr regions. Also triple junctions can occur. 

Therefore, we consider a decomposition of A4 in time-dependent regions ..., 

t G [0, T], separated by curves Ti(t),..., Tjsf c (t). Each curve is parameterized by a time-dependent 
function Xi (. ,£) : If -G M 3 , where I{ is a one-dimensional reference manifold for i = 1,... ,Nq- 
Similar to the case of one curve, we set V<$>^ — and i7 m, i = (xi) s x The geodesic curvature 

AM,i an( i lh e normal curvature are given by = (%i)ss -vM,i an( i = (%i)ss • All 

quantities are time-dependent. 

We define a piecewise constant image approximation by u(.,t) = Ylk=i c k(t)Xn k (t)i where 
Xn k (t) Is th e characteristic function on the set £!&(£) and the coefficients c^(t) are computed by 


Cfe(i) = 


Aft) 

Aft) 1 ^ 


(17) 
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i.e. they are set to the mean of uq in 

Let Xi (., 0), i — 1,..., Nc, be parameterizations of given curves 1^(0) C Ad. We have to solve 
the following scheme for t G (0, T\: Find Xi( ., t) : Ii —> M 3 , ., £), ft$^(., t) : —>► M such that 


0?i)* • + Fil (18a) 

(xi)t • = 0, (18b) 

(^z)ss = “F (18c) 

hold for i = 1,..., A/^. The external force F{ is defined for f G M by 

Fj(f) = fi + A fe +(j)(it 0 (^) - c fc +(j)) 2 - A fc -(j)(« 0 (5) - c fc -(j)) 2 , (19) 


where k + (i), k (i) G {1,..., Nr} are indices of two regions, such that i7 m, i points from Yl k -(i) 1° 

^4+(i) • 

In the experiments, presented in Section |4j we always consider the case fi — 0 and X k = A for 
all k = 1,..., Nr, i.e. all segmentations in our demonstrations can be performed with only two 
weighting parameters a and A. The external forcing term is 

Fi(x) = A[(uo(®) - c k+{i} ) 2 - (uo(x) - c k - (i} ) 2 \. (20) 

We also allow open curves, i.e. curves with dTi(t) ^ 0. Since we consider interface curves, i.e. 
each curve Ti{t) separates two different regions an( i ^k~(i)i we exclude free endpoints. This 

means, we exclude the case that a curve ends in Ad without meeting another curve at its endpoint. 
Further, we consider only smooth, compact surfaces Ad without boundary. Thus, endpoints of 
open curves are therefore part of triple junctions denoted with A& G Ad, k = 1,, Nr- 

For each k G {1,..., A^r} let the integers 4 ,i>4,2>4,3 £ {1 ,...,Nq} denote the indices of 
curves Y ik V l = 1,2,3, 4,l ^ 4,2 ^ 4,3 ^ 4,1 with parameterizations x ik l : I ik l = [0,1] -> Ad, 
such that 

x ik,i(Pk,l) = x h,2 (Pk,2) — x ik,3 {Pk,3^ ~ 4 k , 

where G {0,1} corresponds to the start or endpoint of the curve 4,0 Z = 1, 2, 3. 

At the triple junctions A k , k = 1 ,..., Nr, an attachment condition and Young’s law need to 
hold: 


the triple junction A k does not pull apart, 

(21a) 

3 

E(- 1 ) Pfc ’ i = °. 

1=1 

(21b) 


where T{ kl := (xi kl ) 8 is a tangent vector field at Ti kl C Ad, l = 1,2,3. We refer to |2[ for the 
planar case of evolution of curves with triple junctions. 


2.4 Vector-valued Images 

The segmentation method can be easily extended to vector-valued images such as color images. 
Only the external force F need to be adapted compared to scalar images. The adaptations of F 
can be done similarly as for planar images [8], (4j. 
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Let uq = (^ 0,17 ^ 0 , 2 ? ^ 0 , 3 ) : A4 M 3 represent a color image in the RGB color space. The three 
components of the vector-valued image function represent the red, green and blue color channel. 
The external force Fi, i = 1,..., Nq, in (20) has to be modified. We therefore set: 


3 

Fi ~ [( u 0,j — c k+(i)j) ~ (' U QJ ~ c k~ (i),j ) ] 5 

3 = 1 


where A j weights the j-th color component, j = 1,2,3. The coefficients c^ 
k = 1 ,..., Nr, are given by setting to the mean of uqj in 

Another color space, we will use for segmentation of color images, is the chromaticity-bright ness 
space. Let vo = uo/\\uo\\ be the chromaticity and 60 = \\^o\\ the brightness of an image with image 
function uq. We modify Fi to 

Fi = A c [||V 0 - 4+(i)H 2 - l|V 0 - 4-(i)H 2 ] + A b [(60 - h+(i)) 2 ~ (b 0 - fyfc-(i)) 2 ] , 

where A c weights the chromaticity term and Xr weights the brightness term. Here, is a 
normalized mean of vq in the region (see |4] for details) and b & is the mean of in 


2.5 Restoration of Images on Surfaces with Edge Enhancement 

Searching for a minimizer of the Mumford-Shah functional for images on surfaces Q involves both 
a set of curves T and an image approximation u. The segmentation technique presented above uses 
piecewise constant image approximations to divide an image into characteristic regions of similar 
image intensity or color. For image restoration, more details of the original image uq should be 
kept; the piecewise constant approximation would be a too large simplification. Therefore, we aim 
at approximating uq : A4 R by a piecewise smooth function u : A4 R. 

We propose to first perform a segmentation of the image with piecewise constant approxima¬ 


tions by solving the evolution equations (18) with (21) in case of triple junctions. This is followed 


by a restoration using the already identified regions. By denoising the image in this way during a 
postprocessing step, the edges in the image will not be smoothed out if the curves T match with 
these edges. 

We thus consider piecewise smooth approximations u\Q k = u where : Qk [0,1] is a 
smooth function defined on the region = Nr. For smoothing the image in the 

regions we derive surface partial differential equations from the Mumford-Shah functional. 
Since the curve set V has already been determined, we can fix V in the functional 0 and consider 





variations of u of the form u + ev, for v : M ^ W and e > 0. We compute 


d 

de 


e=0 


E Mb (u + ev, T) = lim ~(E Mb (u + ev, T) - E Mb (u, T)) 


d A 


= lim - ( f {2e\7 mu ■'y mv + <?\\^ mv\\ 2 ) 

e_y0 e \JM\T 

+\ / (2e(u — uo) v + e 2 v 2 ) dT 

Jm J 


= 2 / Vm u ■ ^m v d^4 + 2A / (u-uo)vdA 
Jm\ r Jm 

Nr . 

= 2V / {y M UkMV + K u k ~ Uo) v) dA. 

k=l ^ k 


For a stationary solution, we search for a function u satisfying 0 = ^ | f=(| U MS (u + ev, T). This 
leads to 


0= / (V MUkMV + K u k - uq) v) dA 
Jn k 

for each k — 1,, Nr and an arbitrary function v which is smooth on fl/.. Using an integration 
by parts formula [13], we obtain 


0 = / (~A M Uk +\(u k -u 0 )) vdA+ / V MUk-iJkvds 

J Ofc JdQk 


( 22 ) 


where A m is the Laplace-Beltrami operator and Pk(p) is a unit outer normal vector on dfi& in 
T^M, i.e. it is tangent to the surface for each p G dfife C M but normal to dfi& in P- Since M is a 
smooth, compact surface without boundary, the boundary dfi& of the region fi& consists of one or 
more curves T^, i G {1,... ,Nq}- Thus, locally, pk is ±z/Since v is arbitrary chosen, we have 
to solve the following surface partial differential equation with Neumann boundary condition for 
k = 1 , • • •, Nr: 


~^^M u k + u k =u 0i 

V M u k • /4 =0, 


in fife, (23a) 

ondfife. (23b) 


The smoothing effect is due to the Laplace-Beltrami operator. We can control the smoothing 
extent using the weighting parameter A > 0. The smaller A, the larger is the denoising. The 
larger A, the closer is the approximation to the original image. The Neumann boundary condition 
provides that edges in the image are not smoothed out. 

Vector-valued images can be smoothed by considering each component individually like a scalar 
image. 


3 Numerical Approximation 

3.1 Finite Element Approximation of the Image Segmentation Scheme 


We introduce a finite element approximation for the scheme 
junctions. The evolution equations (18a), for i = 1, ...,iV’c', can be interpreted as a weighted 


with © in the case of triple 
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geodesic curvature flow with external forcing terms. Therefore, we make use of the finite element 
scheme for geodesic curvature flow developed in |3j for closed curves and generalize the approach 
for possible open curves with triple junctions and for image segmentation problems. 

In order to formulate a finite element scheme, we first introduce a spatial and time discretiza¬ 
tion, introduce discrete function spaces and discrete inner products. 

For i — 1 , ...,iVc', let 0 = q z 0 < q\ < ... < q l N — 1 be a decomposition of the interval 
li — I — [0,1]. If Ti is a closed curve, we make use of the periodicity Ni = 0, Ni + 1 = 1, 
—l = Ni~ 1, etc. 

We introduce the following discrete function spaces 


W h : = 


■= {Oi, • • •, Wc) e [Ci 1 ! M )] JVc : is linear, Vi = 1 ,, N c , j = 1,..., v} , 


(24a) 


V h ■= € [C^M 3 )]^ 0 : W kil (p kt i) = ffi kt2 (p k ,2) = fji kia (p k #),Vk = l,...,NT, 


Vi\[ q i_ 1>q i] is linear, Vi = 1 ,..., N c , j = 1,.. •, v} • 


(24b) 


The attachment conditions for triple junctions are incorporated in the definition of the space 
V h . 

A basis of W h is given by functions Xi,j ■= l, • • •, ( Xi,j)N c ) € W h , where (; Xi,j)k(qf ) := 

for i, k = 1 ,..., iVc, j = Jo, • • •, ATj, l = j$ ,..., iV fe , where j l 0 = 1 if If is closed and jjj = 0 

else. 

Further, let 0 = to < t\ < ... < tM — T be a partitioning of the time interval [0, T] into possibly 
variable time steps r m := t m +1 — t m , m = 0,..., M — 1. Let X m (Xf 1 ,..., X™ c ) G V 71 be an 
approximation of x(.,t m ) = (xi(., t m ),..., xw c (.,t m )). Let r m = (r^,..., T 7 ^) denote the 
image of X m . In case of triple junctions let jk : i G {0, denote the index of the corresponding 

curve endpoint such that qk = 1 ,..., iVj 1 , Z = 1, 2, 3. 

For scalar or vector functions u = (izi,..., un c ), v = {v\ G [L 2 (/,R( 3 ))] iVc , the 

L 2 -inner product (•, -) m over the current polygonal curve network T m is given by 


r r 

{u,v) m := / u.vds := V" / 
^r m i=i Jl i 


Ui-Vi ll(X z m )p|| dp. 


(25) 


We follow the ideas of |2 , |3j, and define a mass lumped inner product (•, *)^ for piecewise 
continuous functions i/, = (ixi,..., an( i ^ = Cfl,..., vn c ) by 


1 N C ^ 

mi [(“<■”<) (i«r)+ 

*=1 j = l 


(26) 


where := hm e ^ 0 ,e>o^(< 7 j ± e) and h™ i := ||Xf l (gp - Xf l (^-_ 1 )|| > 0 is the distance 

between two neighbor nodes. 

Let h m := max^=i v ..,N c ,j=i,...,Ni h m _i denote the maximum distance between two neighbor 

Jrf~2 

nodes of the polygonal curves. Let X m G V_ h be a given parameterization of the polygonal curve 
network T m satisfying the following assumption: 
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(A) The distance between two neighbor nodes of T m is positive, i.e. h m ._ i > 0 for i = 

2 


1 ,...,N C , j = 1 and X^(q l j+1 ) ± for i = 1,...,N C and j = if 

dF™ = 0 and j = 1 ,..., W - 1 if dT? ± 0 . 


We define dig 1 = 


to, 


<S>,N C 


), by 


WM) = **(*?$)) = 


V*(Xr(q})) 

llV$(Xr(q}))ll’ 


for i = 1,..., Nc and j = j'q, ..., A^, i.e. uJJh approximates z?<^ at time t m . Further, the tangent 
vector field is approximated by cJ™ = We set 


if rf 1 is closed, or if T™ is an open curve and j 7 ^ 0, A^. For closed curves, we make use of the 
periodicity N{ = 0, Ni + 1 = 1 and — 1 = Ni — 1. For the endpoints of an open curve, we define 



xr(Ql)-x?(qb) 

\\X™(q\) - X™(qti\l 




Xliq^-XTiq^) 


Furthermore, we define u 7 ]^ — (k% ?1 ,..., Nc ) by 




Thus, (2 7 J^ i approximates UM,i at time t m . 

The assumption (A) is necessary, such that is well-defined, see also |3|. 

We define a discrete analogue to the space V ,^ by 

Vl = {ff£V h : rJi.uj% A = 0, i = 1,..., N c ] . (27) 

We now propose the following discrete scheme: Let X° G V h be a given parameterization of 
a polygonal curve network T 0 . We assume that the initial nodes Xf(qj) lie on the surface M. 
Further, we assume that assumption (A) holds for X m , m = 0,..., M — 1 . 

For m = 0,..., M - 1, find 5X m+1 G V% and G W h such that 

ftxm +1 

(- ,x^) h m - °(^\x) h m =(F m ,x)L Vx e w h , (28a) 

An 

( K Xi +1 ^MiV)m + (V s <5X m+1 , V s ? 7 ) m = — (V s X m , W s ff)m, Vr/GZ|, (28b) 

where F m = G . with F™, i = 1..... Nq . is the piecewise linear function 

uniquely given by 

*T(?i) = A[(«o(j?T(9})) - XX - (?;)) - cT-aX ( 29 ) 

where are approximations of the coefficients c*.±(^ at t m , cf. We will later state how 

the coefficients c™, k = 1 ,..., Afe, can be computed. 
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Having found 5X m+1 E V%, we set X m+1 := SX m+1 + X m E V h . 


Before we proceed to prove existence and uniqueness of a solution of the scheme (28), we state 
some very mild assumptions. 


(.4.1) Let * € {1,..., Nc}. If dT^ 1 = 0, we assume that dimspanju;^, ti (q l j),^ t i(qj)}fh = 3. 

(^ 2 ) For each k E {1, N T }, we assume that dimspan{{a% ( {q] k ' l )}j={ 1 }f =1 = 

3. 


Theorem 1 . Let the assumptions (.4), (.4i) and (A 2 ) hold. Then there exists a unique solution 
(5X m+1 ,K^ +1 ) € V_$ x W h to the system 


Proof. The system (28) is linear. Therefore, existence of a solution follows from its uniqueness. 


To prove uniqueness, consider the following system: Find {X, km} £ Y% x W h such that 

(X,X^M)i - = 0, Vx € w\ 

(KM*M,rj)m + <V*X,Vs*l) m = 0, VrfEVl 

We obtain choosing % = K M £ W h in ( |30a ) and rf = X E V% in (30b) 

(ft-M i K'Xt)m ~F (^s-^C^s-^Om = 0. 

From this equation, we conclude km = 0 and X = X c for a constant X c = (Xf,. 

(R 3 ) Nc with X? i = X? 2 = Xf k 3 for all k E {1,... , N T }. Further, X? E M 3 satisfies 


(30a) 

(30b) 


X Nc) £ 


(31) 


for all i = 1,..., Nc and j — jg,..., X^, since X G V%. Inserting n = 0 and X = X c , (30a) reduces 
to 

(X c ,x^)^ = 0, Vx G w h . (32) 

We now choose x = Xi,j ^ W h , i G {1,..., X^}, j G {j'q, ..., Ni} in (32). This yields 


xf.uZMj) = o. 


We conclude Xf = 0 using (31), (33) and the assumptions (A\) and (.A 2 ). 


(33) 


□ 


3.2 Solution of the Discrete System 

Let N — X*, with X* = X^ for closed curves and X* = Ni + 1 for open curves. We make 

use of a small abuse of notation and consider functions in W h as elements in R N and functions in 
V r _ h as elements in 

^ — {(^ 1 ? • • • 5 Z Nc) ^ ) • [ Z ik,l\jk,l — [ Z ik, 2 \jk ,2 ~ \ Z ik,3\jk,3l ^ — 1 : . . . } , 

where Z{ G (M 3 )^ and [zf[j G M 3 is the j-th component of the vector Z{. Functions in V% are 
considered as elements in 

= {(zi,.. .,z Nc ) G X : [zi\j = 0, i = 1,..., N c , j = ... ,X*} , 
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with j'q = 0 for open curves and j'q = 1 for closed curves. Let P$> : (R 3 )^ -+ X$ denotes the 
orthogonal projection onto the space X$. 

In order to state a matrix formulation for the discrete system (28), we introduce the following 
matrices 



/ M 1 •• 

0 \ 


("l, ■■ 

■ 0 \ 


(A 1 •• 

• 0 \ 

M := 



, Nm '■= 



, A:= 




( 0 . 

.. M n c ) 


0 . 

• • J 


V 0 .. 

.. A N ° ) 


where M l G R x N * ? G (R 3 )^ x ^ and A 1 G (R 3x3 )^* xA/ **, i — 1,..., Afc, are defined by 

Mji •— (XiJiXi,l)rrv M.)jl := (XiJi Xi,l ^M.)rm ^jl := (^sXiji ^sXi,l)m M3, 
where M 3 denotes the 3x3 identity matrix. We define b m = ( 6 ™, •••, ^ by 


b T = , +JV+ with b i?j '■= +T> Xij)m> 1 = r • • • , N C , j = Jo, • • • , Ni. 


(34) 


The discrete system ([28]) can be rewritten into the following matrix-vector formulation: Find 


K m+1 € M 7V and 8 xm+l € x s , such that 


-0T m M Nj^Pq 

Pc£> Nj\4 P§> AP$ 


771+1 

k m 

5X m+l 


r b m 

1 m u 

—P$AX r 


(35) 


holds, on assuming that X° £ X. 

Since M is non-singular, we can apply a Schur complement approach and obtain 


ro+l _ 
K M ~ 


(pqAPq + -XPqNmM-'nLPq ) SX m+1 = -P^N M M~ l b m - P<pAX m . 


CTTrr 


—M~ l (n t m Ps, 8X m+1 - T m b m ) , (36a) 

Tm ' ' 

(36b) 


Since P<$> is a projection to a subspace of (R 3 )^, the system matrix of the linear equation (36b) 
is singular as a mapping of (R 3 )^ —> (R 3 )^. However, considered as a mapping of X$ -+ X$ it is 
non-singular if the assumptions (A\) and (A 2 ) hold. 

Since the system matrix is sparse, ( |36b| ) can be efficiently solved with linear effort using an 
iterative solver (with possible preconditioning) or using a direct solver for sparse matrices. In the 
examples presented later in Section [4j we use the UMFPACK algorithm 1121 (direct solver) as 
MATLAB built-in routine for sparse linear systems. 

In the following, we will use the abbreviation X™ := 


3.3 Topology Changes 

Our scheme is based on a parametrization of evolving curves on surfaces. Topology changes con¬ 
cerning the curves are not automatically handled which is often considered as the main drawback 
of parametric methods. 

There are different topology changes that can occur: A curve can split in two sub-curves 
(splitting), two curves separating the same regions can merge to one single curve (merging), two 
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curves separating different regions can touch and a new curve and two new triple junctions occur 
(creation of triple junctions) and a curve can shrink and has to be deleted. 

The latter can be simply detected by considering the length of the curve. If the length of a 
curve is smaller than a predefined tolerance, it will be deleted. The other topology changes will 
occur if two points from different curves or two points from one curve which are not neighbors 
have a small distance. A simple comparison of all nodes would lead to an effort of 0(N 2 ), where 
N is the total number of nodes of the polygonal curves T™,..., . Since we can solve the 


linear equation of our main algorithm (36b) efficiently with linear effort, a sub-algorithm to detect 
topology changes should not result in a too large computational effort. 

For curves in the plane, we extended an efficient method to detect topology changes [4] which 
was originally developed by Mikula and Urban |2l], see also |1|. The method to detect topology 
changes is based on an artificial background grid covering the image domain and consisting of a 
finite set of arrays (squares). If two nodes from different curves or two nodes from different parts 
of one curve lie in the same array of the virtual background grid, a topology change likely occurs 
close to the two points. Using this method, the effort to detect topology changes is O(N). 

In principle, the idea of an artificial background grid can be extended to detect topology 
changes involving curves on surfaces. One can construct a 3D background grid around the surface 
M. Again, we can check whether two points of different curves or different parts of one curve 
belong to the same array of the background grid. 

However, using the Euclidean distance in R 3 to detect topology changes can lead to false 
detections: Surfaces exist where two points can have a small Euclidean distance but their geodesic 
distance , i.e. the length of the smallest curve on the surface connecting the two points, is large. 
In such situations, a topology change does not occur. The geodesic distance would be a better 
indicator for detection of topology changes compared to the Euclidean distance. However, the 
computation of the geodesic distance between two points on a surface is very expensive from a 
computational point of view and cannot be used in a sub-algorithm in practice. 

Let a > 0 be the grid size of the cubes of the 3D background grid. For extending the method 
to detect topology changes from the planar case [4] to the case of curves on surfaces, we have to 
choose the grid size a small enough to exclude such wrong detections as described above. 

For p G M let TpM denote the tangent space and NpM — {TpM) 1 - the normal space. Let 
NM = {(p, n) : p G M, n G NpM } denote the normal bundle. 

For the smooth, embedded hypersurface, we consider the map 


E : NM 


(p, n) p + n. 


Theorem 2 (Tubular neighborhood theorem). Every embedded hypersurface M of R 3 has a tubu¬ 
lar neighborhood, i.e. a neighborhood U C R 3 that is the diffeomorphic image under E : NM —> R 3 
of an open subset V C NJ W of the form 

V = {(p, n) G NM : ||n|| < 6(p)} , (37) 

for some positive continuous function S : M —» R. 

Proof. See [201, Chapter 6, Embedding and Approximation Theorems, or |[l9|, Chapter 4, Vector 
Fields and Differential Equations. □ 


For images on surfaces, we assume the surface M to be a compact, embedded hypersurface. 
As a consequence, set 5o = min {5(p) : p G M} > 0. For each p G M the intersection Bs 0 (p) D M 
is simply connected, which is a consequence of the fact that E\y : V —>► U is a diffeomorphism. 
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Thus, the key idea when extending the algorithm from the planar case to curves on surfaces 
is to choose the grid size a of the auxiliary 3D background grid small enough with respect to Sq, 
such that points from two different parts of the surface (with nearly opposite normal vector 
cannot lie in one array of the grid. 

Topology changes are now detected as follows: 

• Construct an underlying 3D grid with grid size a with ay/ 3 < Sq. Note that the intersection 
of a grid element (=cube of grid length a) with Xi is simply connected. 

• Mark the grid elements with the indices of the curves and the mesh points: We successively 
consider the mesh points X™, i — 1,..., Nq , j = Jo, • • •, A^. If the corresponding grid array, 
in which X™j lies, is empty, the grid is marked with (i, j). 

• If a grid array is already marked with (ii, ji) and if Xf 1 - and X™ ^ are no neighbor nodes, a 
topology change is detected. 

• Since X™ an( i X™j 1 ma y n °l be the P a i r with the smallest distance, we consider a few 
neighbor nodes around X™ an( i X™^ • L e t X™ and X™^ be the pair with the smallest 
Euclidean distance in these two small groups of nodes. They can be found by a pairwise 
comparison, which is not computationally expensive since only a few nodes are involved. 

The topology changes splitting, merging and creation of triple junctions (see explanations 
above) are distinguished as follows: 

• If i = ii, a splitting of the curve T™ is detected. 

• If i / ii, we consider the regions separated by T™ and T™: If k + (i) = k + {ii) A k~{i ) = k~{i\) 
, or alternatively k + (i ) = k~{i\) A k~(i) = k + {i\) holds, a merging occurs. 

• Otherwise, a creation of a new contour and a creation of two new triple junctions occur. 

After having detected and identified the topology change, the curves need to be adapted near 
X™ and X ™^. This involves changing the neighbor relations, changing curve indices in case 
of merging or splitting, and creation of a small new contour with a few nodes in case of triple 
junctions. Details are given in j4j. 

In case of triple junctions, a new curve is created. When creating new nodes, one has to ensure 
that these nodes lie on the surface Xi. In the next section, we describe how nodes can be efficiently 
projected to the surface. 

3.4 Additional Computational Aspects 

Triangulated surfaces In practical applications, a smooth function <f> : M 3 —>► R, such that M 
is the zero level set of $, is usually not provided. Moreover, a surface M is typically given as a 
triangulated surface instead of a smooth surface. 

Therefore, we assume that Xi is a union of triangles of a triangulation T h , i.e. Xi = IJcr h eT h (jh ' 
Note, that the function <E> was only needed to compute Normal vectors to the surface can now 
be easily computed for each triangle a h . For a point p on a curve T m C M, we first need to assign 
p to a triangle a h G T h in which the node lies, to compute n$(p). Further, the color data uq is 
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often piecewise constant and uniquely given by its value on the triangles. To evaluate uq (p), we 
also need to assign the node to its corresponding triangle. 

For each simplex, we can project a vector in R 3 to the simplex plane and can use barycentric 
coordinates to determine if the projected node lies inside the triangle. Surfaces are often composed 
of 10 5 to 10 6 triangles. Therefore, for a given point, finding the corresponding triangle in which 
the point lies results in a high computational effort if no additional knowledge is used. 

For mn — 0 and a curve T™, i G {1,..., Nq}, with nodes X™-, j = j'q, ..., A^, we perform 
a global search only for X™^. For j > j’q, we consider first the simplex to which X^-_ 1 has 

been assigned. If the node X™ is not located in the same simplex, we start a search considering 
successively the neighbor simplices. For m > 0, we can assume that a node has moved only slightly 
on the surface from step m — 1 to m. Therefore, we start the search using the triangle to which 
the node was assigned in time step m — 1. Consequently, a global search has to be performed only 
Nq times at the beginning of the segmentation. 

After the linear system ( |36b| ) has been solved, some of the nodes may not lie exactly on the 
surface. For smooth surfaces (like spheres, tori, etc.), the nodes stay very close to the surface 
if small time steps are used, see [§]. However, for triangulated surfaces, a reprojection onto the 
surface is necessary since V<$> is not continuous. A reprojection onto the surface does not result in 
an additional computational effort: In the next time step, we need to compute = n<$> o x and 
need to evaluate uo again for each node. For both, we have to determine again the closest triangle 
for a point. As described above, this is done by projection of the original node to the triangle 
plane and by using barycentric coordinates. I.e. we already need to determine a projection of the 
original point to the corresponding triangle. 

Computation of regions and coefficients For the external forcing term, we need to determine 
the regions approximations of k = 1,..., Nr. The regions £1™ are separated and thus 

determined by the union of discrete curves T m — T™ U ... UT^. Further, we need to compute the 
coefficients c™ which are the average color values of the image function uq in the corresponding 
regions. 

For m — 0, we need to assign each simplex a h G T h to a region fi 3 . For m > 0, we need to 
update the assignment only in a neighborhood of the curves. 

Let j = 1,2,3, denote the vertices of a triangle a\ l . We assign the simplex to a region 

if its center p a h = (p a hi + p a h 2 + P c r h ^)/^ belongs to the region. In the rare case, that p a h 
lies directly on a curve T™, it is assigned either to or For image segmentation, we 

do not apply any special treatment to simplices which are truncated by a curve. 

For a simplex a h close to a curve with center p a h , we can search for the closest node X™j an d 
consider the sign of (p a h - X™) 

For m — 0, we also need to consider simplices which are not close to a curve. The direction 
cannot be used for remote simplices if the surface Ai is curved. Having assigned a small 
band of simplices around the curves, the remaining simplices can inherit the region index by using 
the neighbor relation between the simplices of the triangulation. 

Motivated by these thoughts, we propose the following algorithm for computation of the regions: 

• For all nodes X™, i = 1,..., Nc, j = j’q, • • •, Ni, we consider the triangle cr h to which X™j 
belongs (see determination of the closest triangle described above). If (p a h — X^-) .c3^(^) 


16 





Figure 2: Illustration how triangles are assigned to a region. 1st sub-figure: image and initial 
curve. 2nd sub-figure: small band after no = 4 steps. 3rd sub-figure: regions colored with mean 
brightness value after assignment of all triangles. The surface data is from the Stanford Computer 
Graphics Laboratory, cf. 
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is positive, the simplex is assigned to otherwise to The indices of neighbor 

simplices of a h are stored in a list. 

• We consider the simplices of the auxiliary list, which have not been assigned to a region yet. 
For a simplex <j h of the list, we search for the closest node point X™ and determine a region 
index using the sign of ( p a h — X ™). ujjfo(qj). We store the neighbor simplices of a h in a new 
list. 


• Having assigned all simplices of the current list to a region, the current list is deleted and 
the simplices of the new list are considered. We repeat the procedure no times. Afterwards, 
a small band of simplices around each curve is assigned to regions. 

• For m — 0, the remaining simplices are considered successively by using again lists of neighbor 
simplices. After the step no, we do not determine the closest node X™-. A new simplex 
inherits the region index directly from its neighbor. 


Figure [2] illustrates the assignment of regions for simplices of a triangulated surface. It shows 
the Stanford Bunn^from the Stanford Computer Graphics Laboratory, cf. 130J, and an image with 
three small discs on its surface. We used no = 4 levels of neighbor simplices to assign the simplices 
of a small band around the initial curve to one of the two regions separated by the initial curve. 
The remaining simplices (marked with dark color in the second subfigure) are finally assigned to 
a region by heritage of the region index. 

The final coefficients are computed as follows: Let C™ u ^\a h denote the sum of the 


color data and n™ the number of simplices belonging to The coefficient for the Chan-Vese 
model is computed by setting c™ — C™/n™. For color spaces like the CB space, we need to make 
use of normalized means for some components of the color j4j. 

For mn > 0, we need to update the coefficients only close the curve, i.e. we consider only the 
simplices of a small band around the curves (using again no levels of neighbor simplices around 
the curves). If a simplex a h changes its region assignment from fi] 71 to fi™, we set 


rtf = rtf + 1 , 




rvm _ pm 


+ U 0\a h : 


Q7TI _ QTTl 


U 0 1 (jh • 


(38) 


1 https://graphics.stanford.edu/data/3Dscanrep/ 
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Global refinement-coarsening strategy We can perform a global refinement and coarsening 
of the curves by using two thresholds Z max and Z m i n for the average distance between neighboring 
nodes. Let N- 71 be the number of nodes belonging to a curve T™. If |T™| /N™ > Z max , we perform 
a global refinement of the curve by inserting a new node between two neighbor nodes. On the 
contrary, if |T™| /N™ < l m i n , we perform a global coarsening, i.e. we delete every second node of 
the polygonal curve which is not a boundary point. 

When inserting a new node between two nodes X™j an d during a refinement, we first 

compute p — (X™- + X 7r7 +1 )/2 and determine the closest triangle (again the search for the closest 
triangle is done efficiently by starting with the triangle in which e.g. X™ lies). Having found the 
closest triangle a h , p is projected orthogonally to cr h . Consequently, all newly generated nodes lie 
on the surface. 


3.5 Numerical Solution of the Image Restoration Scheme 


The scheme (23) for image restoration can be solved numerically with a finite element approach. 
Again, we consider a polyhedral surface M given by a set T h of triangles. 

The image restoration is performed as postprocessing step using the final regions from the 
time step m = M. The surface thus consists of polyhedral regions := k = 1,..., Nr. Let 
Tfh — {(J h G T h : cr h C denote the set of triangles belonging to and let pkj, j — 1,..., 
denote the vertices of the triangles belonging to 

For each k = 1,..., Nr, we define the following finite element space 


Qtl _ 

O b, . 


:= |>eCW£,R) : 


U | (jh 


is linear \/<j h G 


r k h }. 


(39) 


For piecewise continuous functions u h ,v h : with possible jumps at edges of simplices 

cr h G 7?, we define the mass lumped inner product 


{u\v h ) h :=- 

o 


E i^i Y,(n h .v h )m, 


(40) 


a h er£ 


3= 1 


where \a h \ denotes the area of a h , and as above p a h j, j = 1, 2, 3, denote the vertices of the triangle 

-h 
k 


cr h G Tu and 




lim 

p^p a h d ,p£° n 


u h (p). 


Further, for functions u h ,v h G L 2 (OJf. , R( 3 ) j , we define 

(u h ,v h ) := [ 

Jn 


h - fc d A. 


U . V 




(41) 


We consider the following discrete system for each region k G {1,..., Nf>\: Find u h G such 


that 


^MU h y M v h ) + (u\vT = (u 0 , v y, G si 


(42) 


where A > 0 is a weighting parameter (cf. ©)• 
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Let with 4>ki(Pk,j) = $ij denote the standard basis of S k . Using this standard basis, 

we can identify each element in S k with its coefficient vector in R N k . Further, we define the 
matrices G R N k xN k by 

i,3 = !,■■■,Nit. 

The entries of the matrices and A\ are computed by considering each triangle cr h G T k and 
computing the contribution of a h to the entries corresponding to the indices of its vertices. For 
computing the contribution of cr h to A\ we need to compute surface gradients. 

For that, we consider three nodes Pk,jnPk,j 2 an< d 5 ji, J 2 ,J 3 E {1, . .., A^}, being the vertices 
of a triangle a h in T k - For the ease of notation, we assume j i = 1, J 2 = 2 and js = 3. One tangential 
vector is given by 

Pk ,3 ~ Pk ,2 
1 * \\pk,s-PkM\‘ 

A second tangential vector which is orthogonal to t\ can be obtained by 

Pk, 1 ~ Qk 
\\pk,i-qk\\' 

with q k = p ki 2 + ((Pfc,i - Pk, 2 ) . Tl) Ti. 

We note that <9^0-^ = 0. For d^ 2 (\)\ 1 we consider a curve 7 : [0, ||p^i — q k ||] —>► R 3 , 7 (e) = 
+ e T 2 . The composition cfy x o 7 is given by 


i°#) = 


fe,i “ 9fc| 


The derivative of x in direction T 2 i 


is 


dr^k, 1 = ^0fc,l(7(e))|e=O = 


1 


fe,l Qk | 


The surface gradient is then given by 

n + 1 r 2 = ..f *’ 1 jf fc . |2 . 

Similarly, we compute 2 lo-^ an d 

The discrete equation (42) can be rewritten to the following linear system: Find u h G K lv fc 
such that 

\A h k u h + M k u h = m£u 0 , (43) 


A 


holds, where Uq € M jV (‘ is given by 


(tfo)i = 

with T k j = |cr fe € T k : p k ,j € cr h 
for j = 1 ,..., 


E^£7?. I a/l | (“o)|afc 

KiJ 


Yha h &r, h ^ \ a 


j j — 1) • • • Wfc) 


(44) 


k,j 


|. Using this definition of Uq, we obtain (no, 4>k j) h = (M k U,ojj 
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Note, that the Neumann boundary conditions are automatically incorporated in the scheme 


(43). The finite element approach is based on a weak formulation of (23) which contains the 


Neumann boundary conditions as natural conditions. 

We obtain an element-wise constant image approximation U h by setting 


3 

3 

3 = 1 


U h \ a h = l > 


(45) 


for simplices cr h G T^. 

By solving the diffusion equation on each region independently, the boundaries of the regions are 
not smoothed out. Vector-valued images are denoised by applying the method on each component. 


3.6 Summary of the Image Processing Algorithm 


We propose the following algorithm for image segmentation with postprocessing image restoration: 
Given a set of polygonal curves T° = (T?,..., T° Nc ) and X° = (V?,..., X % c ) with X^U) = T°, 

Xf(qj) G M, i — 1,..., Nc , j = Jo, • • •, A^, perform the following steps for m — 0,1,..., M — 1: 


1. Compute the regions 14™ C A4 and the coefficients c™, k = 1,..., Nr, as described in Section 

fT4l 


Compute b m as defined in (34) by using the coefficients c ™ o f ste p [T| 
X m + SX m+1 by solving the linear equation (36b), see Section 


3.2 


Compute V m+1 


3. Check if topology changes occur, see Section [3~~3{ In case of a topology change, except for a 
pure deletion of a curve, repeat the steps [l] and [2] n su b-times with a step size of r m /n su b and 
execute the topology change when it occurs in a substep. 


4. If necessary, perform global coarsening or refinement as described in Section 3.4 


Having found a final segmentation of the image at time m = M, perform a restoration by 


computing a piecewise smooth approximation of the image function as presented in Section 3.5 


4 Results and Discussion 


4.1 Artificial Test Images 

Images on the Stanford bunny We test the developed algorithm for image segmentation by 
considering artificial images on the Standford bunny (Stanford Computer Graphics Laboratory 
[301). In a first experiment, we consider a gray-scaled image showing three dark discs. This 
example is similar to an experiment presented by Kruger et al. 17 , who use a level set method 


to solve a geodesic active contours model with a balloon force for images on surfaces, whereas we 
use a parametric method for a Chan-Vese like model for images on surfaces. 

Figure [3] shows our image segmentation result using the developed direct, parametric approach 
for image segmentation. We use the parameters a — 2 and A = 50 to weight the curvature and 
external term. Let At = r m denote the time step size. The time step size is set to At = 0.01. This 
example demonstrates topology changes; in detail it shows how one initial closed curve is split up 
in three single curves. The contours at four different time steps and the corresponding piecewise 
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Figure 3: Image segmentation of a gray-scaled image on the Stanford bunny. Demonstration of 
topology changes (splitting). First row: Original image and contours for m = 1,100,140,150. 
Second row: Piecewise constant approximation. The surface is from the Stanford Computer 
Graphics Laboratory, cf. 
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constant approximation are presented in Figure [3j Of course, a level set technique as used in |17 
can handle splitting automatically, whereas we need to detect the change in topology explicitly 
using the method described in Section |3.3| However, our method to detect topology changes is 


efficient, since it has a computational effort of 0(N ), where N is the number of node points of the 
polygonal curves. 

In a second experiment, we demonstrate the creation and handling of triple junctions for a 
curve network on the Stanford bunny, cf. Figure [4| For this experiment, we set a = 1, A = 40 and 


At = 0.01. The possibility of triple junctions is not considered in 17 


Further, we apply the image denoising method described in Section T5 and T5, cf. Figure [5| 


Since the term A M u k in (23) is weighted with 1 /A, the denoised image is close to the original 
image, the larger A is chosen, cf. Figure [5j Setting A = 1000 or A = 10000, the noise is not 
completely smoothed out. Setting A = 0.1 the resulting denoised version is close to a piecewise 
constant image approximation. Setting A = 100 results in a good denoised image. 


Image on a torus We consider an artificial image on a torus to demonstrate that the method 
can be applied on surfaces of arbitrary topology type (for example arbitrary genus). Figure [ 6 ] shows 
a torus with a color image from different viewing angles and the curves at different iteration steps 
during the evolution. The RGB color space is used for the segmentation. As weighting parameter 
for the curvature a — 1 is used, and all three components of the color are weighted equally with 
Ai = A 2 = A 3 = 20. The time step size in this experiment is set to At = 0.0001. 

Two different topology changes occur in this example: Around m = 335, two triple junctions 
and a new curve are created. Shortly after m = 422, one blue curve splits into two single curves. 
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Figure 4: Image segmentation of a gray-scaled image on the Stanford bunny. Curve network 
with triple junctions. First row: Original image and contours for m = 1,42,100,250. Second 
row: Piecewise constant approximation. The surface is from the Stanford Computer Graphics 
Laboratory, cf. 1301. 

4.2 Real Images 

Lip contour segmentation We consider an application where lip contours should be detected 
on given face image data. Figure [7] shows the results where the image segmentation algorithm 
is applied on three sample images of the 3D face scans from the 3D Basel Face Model (BFM) 
published by the Computer Science department of the University of BaseQ see also |26|. The 
images are segmented using the chromaticity-brightness color space with a = 0.25, A^ = 200, 
A# = 20 and At = 0.001. The initial contour is a closed curve placed around the lips. Applying 
our algorithm for image segmentation, we obtain the final lip contours. 

Processing of global Earth observation data Another application is the processing of Earth 
observation data. Global Earth observation data can be interpreted as an image given on a 
sphere. We apply the image segmentation and denoising method on data from the NASA Earth 
Observation data set|^J cf. |24j. We segment an image showing outgoing longwave radiation^} 
see Figure [8]- [lOj The colors represent the amount of outgoing longwave radiation leaving the 
Earth’s atmosphere in one month (here: January 2014). Yellow and orange color represent greater 
heat emission (around 300 to 350 Wm“ 2 3 4 ); purple and blue color represent intermediate emissions 
(around 200 Wm -2 ). 

Figure [8] presents the given image data and the contours at different time steps (rows), observed 
from different viewing angles (columns). For the segmentation we used the RGB color space and 

2 http://faces.cs.unibas.ch/bfm/main.php?nav=l-0&id=basel_face_model 

3 http://neo.sci.gsfc.nasa.gov/ 

4 Imagery by Jesse Allen, NASA Earth Observatory, based on FLASHFlux data. FLASHFlux data are produced 
using CERES observations convolved with MODIS measurements from both the Terra and Aqua satellite. Data 
provided by the FLASHFlux team, NASA Langley Research Center. 
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Figure 5: Image denoising with edge enhancement. Subfigure 1-5 (row-wise): Denoising result 
using A = 0.1,1,100,1000,10000. Subfigure 6: Original image with noise. The surface is from the 
Stanford Computer Graphics Laboratory, cf. |30|. 


set the weighting parameters for the curvature and the external term to a — 1 and Ai = A 2 = A 3 = 
A = 50. As time step size we set At = 0.0005. Several topology changes (splitting and merging) 
occur during the segmentation (cf. e.g. mn — 75). Figure [9] shows the corresponding piecewise 
constant approximations. Figure [10] presents the result of the postprocessing image restoration 
using the parameters A = 100, A = 1000 and A = 10000. 

For demonstration of multiphase image segmentation, we consider a second example image 
from the NASA Earth Observation data set, cf. |24]. We now segment an image showing the 
Earth’s net radiation. The net radiation is defined as the difference between the amount of solar 
energy which enters the Earth system and the amount of heat energy which escapes into space 
during one month (here: March 2014). Red color represents a net radiation around 280 Wm“ 2 , 
yellow color a net radiation around 0 Wm -2 and blue-green color a net radiation of —280 Wm -2 . 

Figure [lT] presents the original image with the contours at different time steps (rows), observed 
from different viewing angles (columns). For the segmentation we used a = 1 and Ai = A 2 = A 3 = 
A = 300 and the RGB color space. As time step size we set At = 0 . 002 . The detected regions 
are not separated by sharp image edges. Here, the detected boundaries are weak edges, i.e. they 
lie at locations in the image where the color smoothly changes from yellow to orange or yellow to 
green, respectively. At m = 71 a merging and at m— 149 a splitting occurs. Figure p~ 2 ] shows the 
corresponding piecewise constant approximation. 

Since the original image has little noise, we add some artificial, Gaussian noise to the image 
and repeat the segmentation using the generated image. As postprocessing step, we smooth the 
noisy image applying our restoration scheme, cf. 


Section 3.5 


Figure [13] shows the image with 
added noise and the denoised image for different parameters A under different viewing angles. We 
compare the results for several choices of A: Using A = 100, the denoised image is very close to 
the piecewise constant image. For A = 1000 and A = 10000, we obtain images with removed noise, 
but which still contain sufficient details of the original data set. The data is not smoothed out too 
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Figure 6: Segmentation of an artificial image showing different objects on a torus. First - sixth 
row: Original image and contours for m = 1,100, 200, 335,422, 600. First - third column: Different 
viewing angles. 
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Figure 7: Lip contour detection with a two-phase segmentation. Initial (first row) and final (second 
row) contours. The surfaces and images are from the 3D Basel Face Model (BFM) of the Computer 
Science department of the University of Basel, cf. |26|. 


strong compared to A = 100. Figure [l4| shows a magnification of a part of the surface. It shows 
the noisy image (left) and the result of the denoising using A = 10000 (right). We observe that 
the noise is well smoothed out. 

We have shown how the developed method can be applied on segmentation of global Earth 
data. The given data need not be a classical image generated by a camera. The data can be any 
data defined on the Earth’s surface, like radiation data as in the examples presented here. 

One may argue that global Earth data can also be processed on a flat 2D domain, with a 
rectangular grid given by a discrete set of longitudes and latitudes. In principle, one can apply 
a two-dimensional image segmentation and restoration method as developed in |4] for 2D images. 
However, performing the image processing using the 2D image has some disadvantages: To map 
the global Earth data to a rectangular 2D image, image boundaries at the poles and at longitude 
0 have to be created. Points which are close to each other on the sphere (like at the two opposite 
sides of longitude 0) are not close to each other in the 2D image. Also topology changes like 
boundary intersection will occur in the 2D image. The coordinates of the boundary nodes on the 
left and the right boundary of the image may not fit, i.e. they lie on longitude 0, but can have 
different latitude values resulting in two different 3D points. Further, the length of curves and the 
area of regions near the poles are differently scaled compared to those near the equator. Polar 
regions always appear larger in 2D images. After a segmentation is performed, the size of the 
segmented regions are often of interest. Therefore, the area of the regions can be easily computed 
in a postprocessing step using the sphere data. This is not directly possible from the 2D image 
due to the different scaling of the polar and equatorial regions. In summary, it is beneficial to 
consider global Earth data directly as images on a surface. 

5 Conclusion 

We presented how images on a surface can be efficiently segmented by curve evolution with a 
parametric approach. Furthermore, we showed how restoration with edge enhancement can be 
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performed as a postprocessing step. We considered extensions of the Mumford-Shah 23 and Chan- 
Vese [91 models to images on surfaces. The velocity of the parameterized curves was restricted to 
the tangent space of the surface, which guarantees that all curves modeled as parametric curves in 
M 3 stay on the surface during the evolution. We introduced an efficient numerical scheme based 
on a method for geodesic curvature flow |3|j. Topology changes can be detected fast with an effort 
of 0(N ), where N is the number of node points of the discretized curves. The applicability of 
the developed schemes on different images and different surfaces has been demonstrated in several 
experiments. 
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Figure 8: Segmentation of longwave radiation data given on the Earth’s surface. First - fifth row: 
Original image and contours for m — 1,20,50,75,150. First - third column: Different viewing 
angles. The original image is from the NASA Earth Observation data set, [24] . 
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Figure 9: Segmentation of longwave radiation data given on the Earth’s surface. First - fifth 
row: Piecewise constant approximation for m = 1,20,50,75,150. First - third column: Different 


viewing angles. The original image is from the NASA Earth Observation data set, 24 
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Figure 10: Denoising of the longwave radiation data 
m — M — 150. First row: A 100. Second row: A 
column: Different viewing angles. The original image 
set, 


24 


using the detected regions from time step 
1000. Third row: A = 10000. First - third 
is from the NASA Earth Observation data 
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Figure 11: Segmentation of net radiation data given on the Earth’s surface. First - fifth row: 
Original image and contours for m = 1,50,71,149,180. First - third column: Different viewing 
angles. The original image is from the NASA Earth Observation data set, [24] . 
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Figure 12: Segmentation of net radiation data given on the Earth’s surface. First - fifth row: 
Piecewise constant approximation for m = 1,50,71,149,180. First - third column: Different 


viewing angles. The original image is from the NASA Earth Observation data set, 24 
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Figure 13: Denoising of the net radiation data (image with added noise). First row: Noise added 
image to be smoothed. Second row: Smoothing result using A = 100. Third row: A = 1000. Forth 
row: A = 10000. First - third column: Different viewing angles. The original image is from the 
NASA Earth Observation data set, [24] . 



Figure 14: Magnification of a part of the image. Left: noisy, original image. Right: smoothed 
version with A = 10000. 
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